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We address the interesting temperature range of a glass forming system where the mechanical 
properties are intermediate between those of a liquid and a solid. We employ an efhcient Monte- 
Carlo method to calculate the elastic moduli, and show that in this range of temperatures the 
moduli are finite for short times and vanish for long times, where 'short' and 'long' depend on 
the temperature. By invoking some exact results from statistical mechanics we offer an alternative 
method to compute shear moduli using Molecular Dynamics simulations, and compare those to the 
Monte-Carlo method. The final conclusion is that these systems are not "viscous fluids" in the 
usual sense, as their actual time-dependence concatenates solid-like materials with varying local 
shear moduli. 
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Traditionally the term "glass-transition" refers to the 
spectacular increase in viscosity when the temperature 
T of a liquid is lowered within a relatively narrow tem- 
perature range. Some authors even proposed to call an 
amorphous medium a "glass" when its viscosity increases 
above 10^^ Pa s For sufficiently high temperatures 
such glass-former are liquids with a finite viscosity. For 
sufficiently low temperatures these systems are solids, al- 
though they fail to exhibit any crystalline order. Here 
we address the range of temperature in between, and 
attempt to clarify the state of matter that is observed 
there. Over the years a qualitative picture seems to have 
emerged: already Goldstein in 1969 proposed that in 
the intermediate temperature range the potential energy 
surface exhibits various minima separated by relatively 
high potential barrier (with respect to ksT). Stillinger 
and Weber expounded this further, discussing the exis- 
tence of "inherent states" . The dynamics is trapped for 
a while in an inherent state and then hops, via a saddle, 
to another "inherent state" Demonstrations of such 
events were presented for example in Q . The aim of this 
Letter is to explore the nature of these "inherent states" 
from the point of view of their mechanical properties. We 
will offer a complementary point of view focusing on the 
elastic moduli of the material. An accepted notion of a 
liquid is a medium that cannot support shear, and thus 
by definition has a zero shear modulus. In this Letter we 
propose that it is advantageous to study the shear mod- 
uli in glass formers, since the notion of viscosity loses 
its usual meaning as a local descriptor of the state of 
flow in this range of temperatures, as is explained below. 
In contradistinction, the shear moduli lend themselves 
to accurate measurement also locally in time throughout 
the temperature range. Moreover, the typical jamming 
that occurs at low temperatures can be overcome when 
computing the shear moduli by employing an efficient 
Monte Carlo method 0] as is described below. Using 
this point of view we will be able to present in this Let- 
ter accurate calculations showing that there exists an in- 
teresting temperature range where the elastic moduli are 
finite when computed for short times, but tending to zero 
when computed for longer and longer times. The appro- 
priate length of time should be quantified with respect to 
the statistics of the (relatively rare) relaxation events in 



which the system jumps from one state to the other. We 
reiterate that viscosity in the local sense does not exist. 

For concreteness and in order to be able to check our 
calculations with respect to others we chose to work here 
on the well studied glass former obtained from a two- 
dimensional binary mixture of discs intcractingvia a soft 
l/r^^ repulsion with a "diameter" ratio of 1.4 [^,0]. The 
particles have the same mass m, but half of the parti- 
cles are 'large' with 'diameter' (T2 = 1.4 and half of the 
particles are 'small' , with 'diameter' ui = 1. The three 
pairwise additive interactions are purely repulsive: 



Mab = e ( — ) , a, 6=1,2 
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where = and cTq^ = (cq + a-i^jl. The cutoff radii 
of the interaction are set at 4.5(Ta(,. The units of mass, 
length, time and temperature are m, ui, r = ax^mjt 
and e/kBj respectively, with ks being Boltzmann's con- 
stant. Ref. [TJ found, using molecular dynamics simu- 
lations in the isothermal-isobaric (NPT) ensemble, that 
for temperature T > 0.5 the system is liquid and for 
lower temperatures dynamical relaxation slows down. A 
precise glass transition had not been identified in (3]. In 

Q it was argued on the basis of statistical mechanics 
that there exists a typical length scale that grows expo- 
nentially fast when the temperature decreases. Associ- 
ated with this fastly growing scale there exists a fastly 
growing relaxation time, such that below a certain tem- 
perature the system is jammed for all practical purposes. 
Here we will shed further light on this phenomenon. 

To measure elastic moduli at any temperature one em- 
ploys the relation between the strain fluctuations and 
the elastic properties of a system. In the frame of the 
isothermal-isobaric ensemble the strain simultaneous cor- 
relation functions are given by: 



(2) 



where Z is the following partition function: 

/AG 
cxp( —)de^xdeyydej:y (3) 

and eij are the components of the strain tensor. AG is the 
difference of Gibbs free energy between the strained and 
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the reference state. In the hnear elastic approximation 
for a 2-dimensional system it takes on the form 



AG* 1 r \ / \2 / / 



)2+4/.ey, (4) 



where A, fi and /i' arc the elastic moduli, and Aq is the 
area of the system in the reference state with respect 
to which the free energy is computed. With this form 
of the free energy Eq. ^ reduces to simp le Gaussian 
integrals yielding the wanted relations [l^l, [lH. Here 
we are interested mostly in the shear modulus fi which is 
defined by: 



1 



4 < ely > Aq 



(5) 



The strain fluctuation in ([5]) can be evaluated using 
numerical simulations. An efficient Monte Carlo method 
aimed at calculating this quantity in NPT ensemble 
was described in In essence, the method introduces 
strains into the simulation box by first defining a square 
box of unit area where the particles are at positions Si. 
Then one defines a linear transformation h, taking the 
particles to positions via = hsi. The area of the 
system becomes the determinant \h\. In order to prevent 
rotations of the simulation box, the matrix h should be 
symmetric. After the application of the transformation 
one performs n standard Monte-Carlo moves, after which 
the transformation h changes according to 
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where Ahmax is the maximum allowed change of a ma- 
trix element and is a random number uniformly dis- 
tributed between and 1. At this point one performs 
additional n standard Monte Carlo steps etc. The strain 
tensor in this procedure is calculated by (l^ : 



e=^ih^'hhh,' -I) 



(7) 



where ho =< h > is the reference box matrix which is 
obtained from averaging h until the time of evaluation, 
and I is the unit matrix. In our calculation we chose 
n = 100 for the number of Monte-Carlo steps between 
strain transformations. Note that the strain e depends 
on time both due to h itself and due to (the possibly 
slowly convergent) time dependent /lo- 
in Fig. [1] we present the results of our calculations for 
the shear modulus /i. We have used 1024 particles in 
the 2-dimcnsional simulation cell. After an initial run of 
5000 Tih sweeps (undisplayed in Fig. [T]), these quantities 
were measured as a function of the number of times n^, 
of h sweeps, for the five temperatures T = 0.1, 0.2. 0.3, 
0.4 and 0.5 and the pressure P = 13.5 (in units of e/af). 
Clearly, for T=0.1 the shear modulus converges to a time- 
independent value, indicating the existence of a real solid. 
For T=0.4 and T=0,5 the shear modulus converges to 
zero, indicating the existence of a fluid. For T=0.2 and 
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FIG. 1: (Color online). The shear modulus computed from 
Monte-Carlo simulations as a function of the number of rih 
sweeps, for different values of the temperature 




FIG. 2: (Color online). Trajectories of measured values of e^y 
in the Monte-Carlo simulations as a function of Uh- Upper 
panel: T = 0.1. Middle panel: T = 0.3. Lower panel: T = 
0.5. 



0.3 the shear modulus did not converge, continuing to 
decrease as a function of Uh, possibly approaching zero 
asymptotically, but maybe stabilizing at a finite value 
after a very long time. 

In Fig. [2] we plot the actual evolution of e^y as a func- 
tion of Tih- For r = 0.1 these values fluctuate around a 
constant value, but for T = 0.3 we note the sharp tran- 
sition between different "constant" values that now fluc- 
tuate in time. The trajectory at T = 0.5 represents the 
almost free fluctuations in the value of exy leading to a 
huge value of (e^y) which is associated with a vanishingly 
small shear modulus. 

To complement these flndings we present now a dif- 
ferent method to compute shear moduli. To this aim 
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FIG. 3: (Color online). Results of Molecular Dynamics Simu- 
lations. Black stars: exact calculation of /x^ . Cyan diamonds 
and green circles: simulation results for ^oo . Magenta squares 
and red triangles: the shear modulus according to Eq. (|12p 



we invoke some exact results from statistical mechanics. 
Two important results are due to Zwanzig [l3l | and to 
Wallace [IJIj relating the shear modulus to other com- 
putable objects. The first object is the so-called "infinite 
frequency shear modulus" fioo which is defined as 



_ ^0 2 

A*oo — J, < cr^^ >, 



Ejy ' 1 (8) 

where the microscopic shear stress is given by: 



+ (9) 



where p.j q, is the a component of the momentum of par- 
ticle i and Fj „ is the a component of the force exerted 
on particle i [l5| . Zwanzig had considered a different 
quantity, which is referred to by the same name " infinite 
frequency shear modulus" but which is only appropri- 
ate for liquids, denoted here /i^. In two dimensions this 
object reads 



1 d , ^d^{r). 



(10) 



For interparticle potentials of the form 1/r" (n=12 in 
([!])) the expression (fTO|) yields: 



n - 2 



(P-pT) 



(11) 



Wallace obtained in 2000 the remarkable result (appar- 
ently uncited even once) : 



(12) 



where the shear modulus fj, is the same as defined by 
(ini). We thus learn that statistical mechanics provides 
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FIG. 4: (Color online). Trajectories from Molecular Dynam- 
ics simulations in the stress-energy plane, for three different 
temperatures. U/N is the potential energy of the system per 
particle. Each point represents an average over a duration 
of 50t, connected with straight lines. The time of simula- 
tion is 1.5 X lO^r. During this time with T=0.2 the system 
fluctuates around a single "solid state". During the same 
simulation time when T=0.3 the system hops between three 
different and distinct "solid states". When T=0.4 the fluctu- 
ations are so large that the shear modulus averages to zero 
during the simulation time. 



us an indirect way to measure the shear modulus as the 
difference between the exactly calculable PT|) (in our case 
with n = 12) and the simulationally accessible /Zoo Eq. 
®- 

We performed Molecular Dynamics simulations in the 
canonical (NVT) ensemble with N = 256 particles in a 
square simulation box. The equations of motion were 
integrated using a third order Gear predictor-corrector 
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algorithm, a constant temperature was preserved using a 
velocity rescaling method [Tsj . at each temperature the 
density was chosen in accordance with the simulations 
results in NPT ensemble as described in 0]. We also 
followed this reference in correcting the shear modulus for 
the frozen-in stress. Our results and the results obtained 
in 0] arc presented in Fig. [3] Wc see that as long as 
the system is liquid (T > 0.4), the Zwanzig calculation 
agrees exactly with our simulation of /ioo , predicting zero 
shear modulus, as expected. For smaller values of T, 
^ao deviates from the liquid value, and the difference is 
the Molecular Dynamics estimate for the shear modulus. 
As before, in the interesting region of temperatures this 
method indicates a finite value of the shear modulus, but 
this again depends on the simulation time, and still needs 
to be interpreted. 

To understand what is the state of matter that is be- 
hind these results we present in Fig. 3] the actual trajec- 
tory of the Molecular dynamics simulations in the stress- 
energy plane for three different temperatures and for the 
same simulation time. We see that at T = 0.2 the sys- 
tem fluctuates around one distinct state, giving us a non- 
zero value of as computed from ([5]). For all that one 
can judge we have a solid with a finite value of its shear 
modulus. Changing the temperature to T = 0.3 we rec- 
ognize that the system fluctuates around one such "solid- 
like" state, but then jumps to another such "solid-like" 
state, and then another such state; for this time of the 
simulation one resolves at T = 0.3 three states each of 
which would have given us a sharply defined shear mod- 
ulus for the time of life of that state. Upon averaging 
over longer periods we get contributions from different 
"solid-like" states and the shear modulus becomes de- 
pendent on the time of averaging, as seen in Figs. [T]and 
m For T = 0.4 the trajectory now fills up an extended re- 
gion in the stress-energy plane, and one can see why the 



shear modulus must average to zero for longer averaging 
times (see also Fig. [5]). Anyway we look at the dynamics 
the emerging picture is the same: between real solid and 
real liquid the system is locally not a "viscous fluid". 
Rather, the trajectory concatenates relatively long pe- 
riod of times where the system behaves like a solid, in- 
terconnected by relatively short period of times where 
the system flows between these states. Clearly, a viscous 
fluid behaves very differently, responding to stress by a 
viscous flow, be it as sluggish as one wishes. Here, most 
of the time, the system does not respond to stress, except 
in the narrow corridors of flow which become rare when 
the temperature goes down and more common when the 
temperature warms up. Of course this does not mean 
that in the sense of long time averaging the notion of 
viscosity cannot be resurrected, but locally in time this 
is impossible. 

The question that remains is how to describe macro- 
scopically this state of matter. In our opinion a promis- 
ing approach was described in [l^ , where all the acces- 
sible stationary points of the potential energy surfaces 
V(ri ■ • ■ rjv) for all the N particles were obtained by solv- 
ing the 37V nonlinear equations dV/dri = 0. Having a 
complete enumeration of the accessible minima, the sad- 
dles and the eigenvalues of the Hessain matrix at the 
saddle (cf. [3 for details) one can hope to analyze the 
relative length of time that the system might stay in each 
of its "solid-like" minima. In that case one can attempt 
to answer the question whether a long time average would 
yield a finite or a zero shear modulus or at any rate de- 
scribe the mechanical behavior of the system for a finite 
time. Such a program is becoming more and more real- 
istic with the increase in computer power, and appears 
unavoidable if one wants to make real progress in under- 
standing the interesting range of temperatures between 
liquid and solid in amorphous glass-formers. 
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